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integration,  the  nested  grid  solution  is  improved  over  the 
coarse  mesh  solution.   Difficulties  arise  at  the  boundaries, 
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I.   INTRODUCTION 

The  effects  of  mountain  ranges  on  synoptic  scale  motions 
have  been  shown  (Reiter,  1971)  to  be  quite  influential.   Yet, 
the  physical  cause  and  mathematical  description  of  these  ef- 
fects are  poorly  understood.   The  large-scale  flow  over  and 
around  mountain  ranges  is  not  properly  observed  or  mathema- 
tically well  modeled.   This  was  quite  evident  to  this  author, 
having  spent  three  years  of  forecasting  weather  for  the  Euro- 
pean-Mediterranean area,  where  the  Alps,  Pyrenees,  and  other 
mountain  ranges  have  a  marked  meteorological  effect.  Due  to 
the  grid  size  of  most  primitive  equation  models,  a  nested, 
finer  resolution  grid  is  required  to  depict  accurately  the 
flow  over  and  around  mountains.   This  is  attempted  in  this 
s  tudy . 

The  general  circulation  model  used  in  this  study  is  a 
modified  version  (Monaco,  1975)  of  the  UCLA  global  predic- 
tion model.   To  reduce  computer  memory  and  run-time  require- 
ments prior  to  nesting  a  finer  resolution  grid,  a  channeling 
of  the  global  model  was  undertaken.   Procedures  had  to  be 
developed  to  take  into  account  the  scheme  C  (Monaco,  1975) 
grid  arrangement  of  the  model  and  the  higher  order  form  of 
the  flux  terms  in  the  primitive  equations . 

In  nesting  a  fine  resolution  scheme  C  grid,  certain  con- 
ditions had  to  be  considered  to  take  into  account  the  stag- 
gered grid  arrangement.   Unlike  the  majority  of  nested  grid 


11 


models  presently  available  which  use  unstaggered  grid  ar- 
rangements, the  indexing  scheme  used  for  each  variable  must 
be  considered  in  determining  interface  boundary  conditions. 

The  solutions  obtained  are  divided  into  three  categories 
a)  coarse  mesh  solution,  b)  uniform  fine  mesh  solution,  and 
c)  nested  solution.   Comparisons  are  made  between  solutions 
with  no  mountain  present,  and  then  with  a  mountain.   It  is 
shown  that  a  number  of  problems  remain  for  the  type  of 
channeling  and  nesting  used  in  the  present  model. 
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II.     MODEL  DESCRIPTION  AND 

INITIALIZATION  PROCEDURES 

The  general  circulation  model  used  in  this  research  has 
been  described  by  Monaco  and  Williams  (1975) .  It  is  a  modi- 
fied version  of  the  UCLA  global  prediction  model  detailed  by 
Arakawa  and  Mintz  (1974) .  Our  model  assumes  the  atmosphere 
is  adiabatic  and  frictionless  and  contains  no  sink  or  source 
terms.  Numerical  integration  is  carried  out  on  a  staggered, 
spherical,  sigma  coordinate  system. 

A.   VERTICAL  COORDINATE 

The  vertical  coordinate  used  in  the  model  is  the  non- 
dimensional  a-coordinate,  fig.  1.   It  is  defined  as  follows: 


P  -  Pt 

o   =  — - — -  (2.1) 


71  ~   Ps  "  Pt 


where  it   is  the  terrain  pressure.   The  pressure  is  given  by 
p  ,  the  surface  pressure  by  p   ,  and  the  constant  tropopause 
pressure  by  p   .   The  vertical  boundaries  of  the  coordinate 
system,  the  earth's  surface  and  the  constant  tropopause 
pressure  level,  follow  from  (2.1),  such  that: 


=  0  at  p  =  p   (=  200  mb) 


o   =   l     at  p  =  ps  . 
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The  boundary  condition  at  the  tropopause  pressure  level, 

dg. 
St- 


and at  the  earth's  surface,  is:   a(=  -rr-)  =  0  . 


B.   HORIZONTAL  DISTRIBUTION  OF  VARIABLES 

The  dispersion  of  local  accumulations  of  wave  energy 
generated  by  gravity-inertia  waves  is  an  important  process 
in  any  atmospheric  numerical  model.   Winninghoff  (1968) 
showed  how  this  geostrophic  adjustment  depended  on  the  dis- 
tribution of  variables  over  the  grid  points.   Five  possible 
distributions  of  the  dependent  variables  <j> ,  T,  v,  and  u  are 
shown  in  fig.  2.   The  original  UCLA  general  circulation 
model,  described  by  Arakawa  and  Mintz  (1974),  used  scheme  B 
for  the  horizontal  distribution  of  variables;  this  research 
model  uses  scheme  C. 

Winninghoff  showed  that  schemes  B  and  C,  in  a  one-dimen- 
sional case,  adequately  simulated  geostrophic  adjustment. 
However,  in  a  two  dimensional  case,  scheme  C  was  shown  to 
be  the  best  lattice  to  simulate  the  geostrophic  adjustment 
process.   In  using  scheme  C,  though,  two  conditions  can  cause 
problems.   The  first  is  due  to  the  two  velocity  components 
not  being  carried  at  the  same  latitude,  causing  some  diffi- 
culty with  Coriolis  force  terms.   The  assumption  of  a  con- 
stant Coriolis  parameter  could  eliminate  this  problem.   The 
second  problem  occurs  if  ^-o/a      is  less  than  or  close  to  one. 
The  quantity   AR  ,  the  Rossby  radius  of  deformation,  is  given 
by 

xR  =  sgs.  ii 
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where  H  is  the  mean  value  of  h  ,  the  depth  of  the  atmos- 
phere.  The  situation  where  ^t>/a   £  1   i-s  rare  and,  there- 
fore, scheme  C  is  used. 

C.   PRIMITIVE  EQUATIONS 

The  primitive  equations  are  written  using  the  orthogonal 
curvilinear  coordinates,   £   and  n  ,  where   £  =  A  (longi- 
tude) and   n  =  <J>  (latitude).   Scaling  factors,   m  and  n  , 
are  defined  as  follows: 


—  =  a  cos  d>   and  —  =  a  . 
m         Y       n 


The  horizontal  velocity  components,  temperature,  and  surface 
pressure,  the  prognostic  variables,  are  governed  by  the  fol- 
lowing primitive  equations: 


3      /   77  v       ,      3      /7TU        x       -      3      /7TV        \       ,      3      /TTG         \  /n     o  \ 

ytfe  u)    +  T^IT  u)    +  ^(^T  u)    +  7a(ET  u)  (2°2) 


.    f       .   /         3         1  3         1\  ,  77r3d)      ,  3lTn  TT       _ 


i_CJL  v)    +  2_(ZEH  v)    +  1_(IX  v)    +  ^-(li  v)  (2   3) 

3tW      ;         3^  n    v;         3n^m    v;    r   3a  4m  v;  ^*J; 

,.  f     ,/      3      1              3      1N,          ,    7Tr3d)     ,    „      3tt.  _     it  _ 

+  [ —  +  (v  TT?-  —  -   u  TT—  — )  ]  ttu  +  — [•r-L  +  c^a   -~— 1  =  — —  F      , 

mn      v      3£   n             3n   nr                m   3n                3n"         mn  n    ' 


Ttfe  cpT)    +   3T(1T  cp°       3TT(-m-  cpT) 


+    (.£-)K  |_ (I£  c    0)  (2.4) 

vp    y       3avmn     p    y 


r  3      ,    77   N       ,      U     377  V3TT.  TTp. 

™a  ^tW    +KIT  +  m^]    +iEQ    ' 
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and 


9tvmn;         8C^n;         3rTm;         9aSmV         u    '  U.3; 


a   =   ^  (2.6) 


6$  =    -7ra5a    .  (2.7) 

A  complete  listing  of  the  symbols  used  in  the  preceding  equa- 
tions may  be  found  in  the  front  of  this  report. 
The  vertical  a-velocity  is  given  by: 

a 
Tra  =  -  f       V  •  (TrV)da  -  a  f^^)  .  (2.8) 

o 
In  determining  geopotentials ,  conservation  of  total  energy 
under  adiabatic  and  frictionless  processes,  and  conservation 
of   9   and   ln6  ,  integrated  over  the  entire  mass  under  adia- 
batic processes,  are  needed.   In  addition,  the  vertical  dif- 
ference scheme  must  maintain  the  property  of  the  vertically 
integrated  horizontal  pressure  gradient  force.   This  is  a- 
chieved  essentially  by  combining  an  alternate  form  of  the 
hydrostatic  equation,  written  in  terms  of  potential  tempera- 
ture (2.9),  with  the  vertical  integral  of  hydrostatic 
equation  (2. 10) : 

3*  =   c     JS5&   H£-)K    ,  (2.9) 

r1  r1 

I     d(<Da)    =    -     /        (iraa    -    $)da    .  (2.10) 

o  o 
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The  system  is  complete  when  the  above  two  diagnostic  equa- 
tions are  included  with  the  primitive  equations. 

D.   ANALYTIC  INITIALIZATION 

The  initial  wind  and  pressure  (geopotential)  fields  are 
defined  mathematically  in  this  model  which  has  some  advan- 
tages over  obtaining  the  fields  from  a  weather  map.   The 
analytic  initial  conditions:   1)  simplify  the  task  of  inter- 
polating and  balancing  initial  fields  from  constant  pressure 
surfaces  to  sigma  surfaces,  2)  minimize  imbalances  between 
wind  and  pressure  fields  thereby  reducing  gravity-inertia 
waves . 

Kaurwitz  (1940),  in  defining  the  initial  velocity  fields, 
solved  the  linearized,  non-divergent  vorticity  equation  for 
the  stream  function.,   The  initial  distribution  of  geopoten- 
tial was  obtained  by  using  this  solution  in  the  forcing  func- 
tion of  the  non-linear  balance  equation  (Phillips,  1959). 

Initialization  of  the  model,  based  on  the  aforementioned 
solutions,  was  achieved  using  the  following  equations: 

2 

\p     =   -Ba     sin  <f> 

$'    =   a2  A(<j>) 

u     =   -   |  |1  =   Ba   cos   <j>  (2.11) 

a    dtp 

v    = I ii  =  o 

a  cos   §    d\  ' 

where  A(<J>)    =  \   (2P.  +  B)B    cos2    <j)    . 
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These  equations  represent  solid  rotation  of  the  atmosphere 
because  of  the  absence  of  vertical  wind  shear.   The  horizon- 
tal or  zonal  wind  distribution  varies  as  the  cosine  of  the 
latitude  and  the  magnitude  of  the  constant,  B.   The  zonal 
wind  distribution  in  the  model  varied  from  u  =  15.2  m-sec" 
at  34  N  to  u  =  13.8  m-sec   at  58  N.   The  initial  tempera- 
ture field  was  determined  in  accordance  to  the  NACA  standard 
atmosphere  as  defined  by  Haltiner  and  Martin  (1957) . 

E.   PROCEDURE  FOR  INTRODUCING  TOPOGRAPHY  INTO  MODEL 

The  scheme  used  to  introduce  mountain  topography  into 
the  model  is  described  in  detail  by  Hayes  (1977).   Initially, 
heights  at  the  surface  were  set  equal  to  zero  everywhere. 
During  the  model  integration,  the  heights  of  the  mountain 
grid  points  were  increased  as  a  function  of  time  until  the 
desired  heights  were  achieved.   The  problem  with  wind  field 
initialization  in  the  vicinity  of  the  mountain  precluded 
starting  with  a  mountain  "already  built".   The  following 
equations  were  used  at  the  mountain  grid  points : 


=  M  sin2  (5jj-)  ,  t  <  t 
P         P 
sfc 

=  M  ,  t  >  t 


(2.12) 


P 
where   <p  ^    is  the  geopotential  height  at  time   t  ,  M   is 

the  final  geopotential  mountain  height,  and   t    is  the 

period  of  mountain  "building".   Values  used  in  this  research 

are : 


*sfc(max)  =  750°  SPm 

t   =12  hours 
P 
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III.   CHANNELING  THE  COARSE  GRID  MODEL 

In  any  numerical  modeling,  considerations  of  computer 
run- times  and  memory  requirements  are  necessary.   This  is 
especially  true  when  developing  and  testing  fine  mesh  grids 
where  the  grid  spacing,  and  therefore  the  time  step,  are 
small.   With  this  in  mind,  an  attempt  to  channel  the  global 
prediction  model  (Arakawa,  1972;  Arakawa  and  Mintz,  1974; 
Monaco  and  Williams,  1975)  in  the  north-south  direction  was 
made  prior  to  nesting  of  a  finer  resolution  grid. 

A0   PROCEDURES  FOR  CHANNELING 

In  the  initial  phase  of  this  research,  the  assumption 
of  no  variation  in  the  north-south  (n)  direction  was  con- 
sidered.  This  led  to  reducing  the  amount  of  north-south 
pressure  grid  points  in  the  model  from  46  to  7.   In  accom- 
plishing this,  the  pole  modification  schemes  for  the 
various  predictive  equations  were  eliminated. 

In  addition  to  removing  the  smoothing  effect  of  the  pole 
modifications,  the  effect  of  another  scheme  which  zonally 
smooths  the  £ -components  of  the  pressure  gradient  force  and 
the  divergence  was  reduced.   This  zonal  smoothing  eliminates 
computational  instabilities  due  to  the  convergence  of 
meridians  at  the  pole.   This  smoothing  is  a  function  of  lati- 
tude, grid  size,  and  wave  number;  therefore,  due  to  the 
channel  being  centered  at  46  N  with  a  width  of  only  3   lati- 
tude, its  effect  was  much  reduced. 
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Although  the  above  routines  might  have  had  a  damping 
influence  on  some  of  the  large  variations  encountered  with 
the  model,  the  elimination  or  reduction  of  these  routines 
did  not  generate  any  oscillations.   Rather,  these  oscilla- 
tions were  generated  by  the  channel  boundary  conditions. 

B.   BOUNDARY  CONDITIONS 

Simple  cyclic  conditions  east-west  were  used,  but  the 
north- south  boundary  conditions  proved  to  be  quite  formid- 
able.  The  main  constraint,  for  the  north-south  direction, 
was  that  of  no  net  mass  flux  in  the  n -direction.   The  two 
main  problems  in  this  regard  were  the  balancing  of  pressure 
gradient  term  to  the  u-velocity  component  and  the  integrated 
mass  flux  in  the  n -direction.   Compounding  the  problem  fur- 
ther were  the  scheme  C  grid  arrangement  and  the  complex  flux 
forms  of  the  advective  terms  in  the  equation  of  motion. 

The  scheme  C  grid  point  did  not  allow  easy  selection  of 
the  boundary  walls,  due  to  the  asymmetry  of  the  variables 
at  each  (I, J)  grid  point.   The  flux  terms  of  the  advective 
equations  involved  a  second  order  Jacobian  scheme  which 
required  specification  of  some  terms  outside  the  walls. 
These  problems  were  implicitly  inherent  in  all  boundary  con- 
siderations . 

For  the  constraint  of  no  mass  flux  in  the  north-south 
direction  to  be  achieved,  the  v-equation  of  motion,  (2-3) , 
integrated  in  the  vertical  and  ^-directions  must  balance. 
The  three  terms  that  must  be  specifically  dealt  with  along 
the  boundaries  are  the  integrated  pressure  gradient,  normal 


20 


component  of  mass  flux,  and  Coriolis  terms.   (Details  are 
given  in  Appendix  A.) 

The  method  used  to  prevent  net  mass  flux  across  the 
wails  was  to  reflect  the  normal  component  of  the  mass  flux 
from  the  immediate  interior  grid  point  by  one  opposite  in 
sign  but  equal  in  magnitude  just  outside  the  wall.   This  was 
accomplished  by  the  following  equations: 

<™¥>I.2+   <™t£>I.1-   °  <3"X> 


(ttv  ^i)T    ,  +    (ttv  -^)t    t         ,,    =   0    .  (3-2) 

v         ml, Jmax         v         m ' 1 , Jmax+1  v        ' 


The  term  (ttv  — ^-)  T  ,  ,-,  is  a  mass  flux  quantity  outside 
the  boundary  wall,  which  was  needed  in  the  advective  terms 
of  the  equations  of  motion. 

Although  this  procedure  assured  no  mass  flux  across  each 
boundary,  it  did  not  assure  zero  net  mass  flux  over  the 
entire  field.   This  is  due  to  the  inequality  between  the 
northern  and  southern  mass  flux  terms,  i.e.: 

(ttv  M.)     +  (ttv  £i)_  _    ,  i  +  0  .  (3-3) 

v    m ' 1 , 1    v    m  y  I ,  Jmax+1  '  v    ' 

This  imbalance  did  contribute  to  total  field  mass  variations 
with  time. 

Along  the  walls  the  balance  between  the  pressure  gradi- 
ent and  the  geostrophic  u-component, 

f  TT  TT  .  d<$)      ■       8TT.  /t  /  \ 

5E  u  =  "  H[TrT  +  aa  TrT]    '  (3"4) 
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was  difficult  to  maintain.   The  gradient  of  terrain  pres- 

sure,   -*—  ,  was  restored  to  its  initial  value  after  each 

time  step.   This  allowed  the  pressure  to.  vary  along  the  walls, 

but  maintained  the  initial  pressure  gradient.   The  geopoten- 

tial  height  gradient  was  controlled  in  a  different  way. 

Since  geopotential  height  fields  are  not  stored  but  calculated 

at  every  time  step,  restoring  the  gradient  ■*—  to  its  initial 

OT) 

value  was  not  straightforward. 

The  method  used  to  specify  y-  was  through  the  parameters 
that   0   was  a  function  of,  i.e.  temperature,  pressure,  and, 
indirectly,  the  velocity  components.   Since  the  boundary  grid 
points  were  not  predictive  points,  the  values  needed  at  the 


boundary  walls  were  provided  by  extrapolation  schemes  using 

3n 


interior  predictive  values.   Since  tt—  =0  at   a  =  1  ,  but 

1         »  A  n  » 


8tt 
increases  with  elevation  as  oa   -~—  decreases,  any  imbalance 

in  the  parameters  used  to  control   $   become  more  significant 

with  height.   With  or  without  a  mountain  in  the  model,  the 

largest  variation  of  the  v-component  always  occurred  at  sigma 

level  1  and  decreased  downward.   This  was  evidently  due  to 

the  large  -^—  term  in  the  v-component  equation. 

With  the  inequality  of  normal  mass  flux  terms  plus  the 

imbalance  in  the  boundary  zonal  flow,  there  was  a  variation 

in  the  total  field  mass,  and  the  oscillations  were  evident 

in  the  surface  pressure  fields.   This  variation  of  pressure 

along  the  north-south  direction  with  time  is  shown  in  Fig.  3. 

The  approximate  period  of  the  oscillations  is  20  hours,  which 

indicates  that  the  forcing  might  originate  from  the  southern 
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boundary  where  the  inertial  period  is  21.5  hours.   This 
might  indicate  a  net  mass  inflow  from  the  south  in  the  early- 
stages  . 

One  attempt  was  made  to  eliminate  these  pressure,  i.e. 
mass,  variations.   The  change  in  terrain  pressure  between 
two  time  steps,  averaged  over  the  entire  grid,  was  subtracted 
from  the  forecast  pressure  fields.   As  shown  in  Fig.  3,  this 
reduced  the  oscillations  substantially  but,  as  shown  in 
Fig.  4,  caused  an  increase  in  the  temperature  variation.  This 
temperature  increase  then  caused  problems  in  the  calculation 
of  geopotential  heights  along  the  boundaries.   The  procedure 
was  discarded,  thereby  allowing  net  mass  changes  over  the 
entire  field  to  occur. 

The  phase  and  period  of  these  pressure  oscillations  were 
dependent  on  the  grid  and  domain  size.   The  oscillations  in 
a  model  with  a  grid  distance  of  4  latitude  versus  1.5   lati- 
tude are  compared  in  Fig.  5.   In  addition,  the  domain  size 
of  the  1.5   latitude  grid  is  one-third  that  of  the  4   latitude 
grid.   The  phase  shift  of  six  hours  between  the  two  grid 
spacing  is  obvious,  but  there  is  also  a  slight  period  change. 
The  smaller  grid  has  a  period  of  approximately  18  hours. 
This  is  the  inertial  period  for  the  southernmost  run  which 
adds  evidence  that  the  mass  variations  are  forced  from  the 
south. 

C.   GRID  ARRANGEMENT  AND  ADVECTIVE  TERMS 

As  mentioned  in  the  preceding  section,  the  scheme  C  grid 
arrangement  and  the  complex  form  of  the  advective  terms 
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compounded  the  problems  at  the  boundary  walls.   The  follow- 
ing is  a  brief  description  of  these  problems. 
1.   Scheme  C  Grid  Arrangement 

In  Figs.  6A  and  6B  the  parameters  enclosed  by  a  tri- 
angle comprise  a  scheme  C  (I, J)  grid  point.   This  indexing 
nomenclature  lends  itself  readily  to  placing  a  boundary  wall 
not  at  one  latitude  but  two.   This  was  used  in  this  research, 
and  proved  to  cause  more  problems  than  it  solved. 

Since  the  boundary  walls  are  not  predictive  points, 
values  of   tt  ,  T  ,  U  ,  V  must  be  provided  to  allow  time 
stepping  of  the  interior  points.   As  discussed  in  section  B, 
the  extrapolated  methods  used  in  this  research  were  not  com- 
pletely satisfactory.   An  example  of  this  is  the  divergence 
calculation  at  Tr-point  (1+1,  Jmax-1)  .   The  l  ^     term  is 

at, 

based  completely  on  predictive  values,  where  the  -^ — - 
term  is  not.   Since  the  same  gradient  of   tt   and  V  is  main- 
tained between  (Jmax-1,  Jmax-2)  and  (Jmax,  Jmax-1),  an 
artificiality  could  develop  with  the  presence  of  a  distur- 
bance along  the  boundary.   As  long  as  an  indexed  (I, J)  grid 
point  is  used  as  a  channel  wall,  this  problem  will  exist. 

An  alternative  channel  wall  boundary  was  considered, 
but  not  tried.   It  was  to  have  had  an  unequal  number  of 
rows  between  the  parameters   U  ,  <j>  ,  T  ,  tt  ,  a   and   V  in 
the  north-south  direction.   This  would  allow  the  U  ,  <{>  , 
T  ,  tt   latitude  row  to  be  predictive  along  both  the  north 
and  the  south  boundaries.   This  would  also  eliminate  the 
problems  of  imbalances  caused  by  the  use  of  extrapolation 
methods  for  providing  values  at  the  walls.   Indexing  does 
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become  much  more  complicated,  though,  and  this  was  the  pri- 
mary reason  it  was  not  attempted. 

2.   Flux  Form  of  the  Advective  Terms 

The  prevention  of  non-linear  instabilities  in  a  non- 
divergent  flow  can  be  achieved  by  conserving  the  mean  square 
vorticity,  and  mean  kinetic  energy  with  time.   The  global 
model  used  in  this  thesis  was  written  to  simulate  slowly 
changing  quasi-geostrophic  motion.   Prevention  of  non-linear 
instability  called  for  a  second  order  finite  difference 
Jacobian.   This  second-order  Jacobian  was  used  by  Arakawa 
(1974)  to  write  the  non-linear  advective  terms  of  the  primi- 
tive equations.   Due  to  its  second-order  nature,  the  mass 
flux  expression  for  a  specific  grid  point  required  advective 
mass  flux  terms  as  far  away  as  two  (I, J)  grid  points.   This 
presents  a  problem  as  the  boundary  is  approached  by  requiring 
the  specification  of  extra  mass  flux  terms  outside  of  the 
channel  wall  (line  Jmax+1   in  Fig.  6A) . 

The  value  of  the  normal  mass  flux  outside  the  channel 
wall  is  a  reflected  value  of  Jraax.   The  mass  flux  at  Jmax, 
though,  is  an  extrapolated  value  and  not  a  predictive  value. 
This  is  due  to  the  scheme  C  grid  arrangement  and  choice  of 
location  of  the  channel  wall.   At  the  southern  channel  wall, 
a  predictive  mass  flux  quantity  at  3=2   is  reflected  out  to 

J-l. 

If  the  channel  wall  is  placed  at  the  latitude  row 
having   $  ,  T  ,  ir  ,  U  ,  and  a    ,  then  a  simplified  flux  scheme 
would  be  required.   Although  no  extra  flux  terms  would  be 

A  ft 

needed  for  predicting  temperature,  another  external   (ttu  — ) 
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term  would  be  needed  at  J=0  and  J=Jmax+l  to  calculate  momen- 
tum flux  centered  on  a  U  point  (Fig.  6B) .  It  would  probably 
be  adequate  to  simply  set  the  following: 


<™t?>i.o-   (™l?>i,l  =  ° 


(7m  All)       .  (7TU  Aik         =  0 
n  y I , Jmax    v    n  '  I , Jmax+1 


To  briefly  summarize,  the  scheme  C  grid  arrangement 
and  indexing  readily  leads  to  certain  boundary  wall  condi- 
tions.  These  proved  to  be  unsatisfactory,  due  to  extrapola- 
tion techniques  for  providing  the  channel  wall  values. 
Another  channel  wall  scheme  was  presented  which  would 
eliminate  extrapolating  values  to  the  wall.   However,  with 
the  integrated  mass  flux  in  the  n   direction  not  being  zero, 
with  the  changing  of  the  flux  scheme  for  the  advective  terms, 
and  with  the  large  indexing  problems,  this  scheme  can  still 
create  difficulties.   Time  did  not  permit  attempting  the 
latter  scheme. 
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IV.   NESTED  MODEL 

A.   SCHEME  C  NESTING 

The  nested  grid  used  in  this  thesis  consisted  of  two 
uniform  grids  with  the  fine  mesh  grid  centered  in  the  chan- 
neled coarse  mesh  grid  (Fig.  7).   Due  to  the  scheme  C  stag- 
gered grid  arrangement,  the  ratio  of  the  coarse  mesh  length 
to  the  fine  mesh  length  is  three  to  one,  which  corresponds 
to  a  reduction  of  latitude  spacing  of  444  km  to  148  km  at 
the  center  latitude.   The  longitude  reduction  is  from  348  km 
to  116  km,  also  at  the  center  latitude.   The  total  number  of 
grid  points  in  the  coarse  mesh  for  one  sigma  level  is  20  x  7, 
while  in  the  fine  mesh  it  is  19  x  7. 

Although  the  fine  mesh  grid  is  placed  at  the  center  of 
the  coarse  mesh,  there  is  a  minimum  number  of  J  rows  required 
to  the  north  and  south  of  the  nested  grid  to  allow  inter- 
action between  the  grids  at  predictive  points.   Due  to  the 
extrapolation  schemes  used  at  the  coarse  mesh  boundaries,  at 
least  two  J  rows  are  needed  to  the  north  and  south  of  the 
nested  grid.   If  the  channel  walls  are  changed  to  consist 
only  of  the   <j>  ,  tt  ,  T  ,  U   latitude  row,  then  this  require- 
ment would  be  unnecessary. 

The  physics  of  the  model  for  the  nested  fine  mesh  are 
identical  to  those  of  the  coarse  mesh  domain.   The  only 
difference  between  grids  was  the  complete  elimination  of  the 
zonal  smoothing  routine  used  in  the  coarse  mesh. 


B.   BOUNDARY  CONDITIONS 

There  are  numerous  approaches  for  specifying  the  inter- 
face boundary  conditions  in  a  nested  grid  model  (Ookochi, 
1972;  Harrison  and  Elsberry,  1972;  Chen,  1973;  Moss  and 
Jones,  1973;  Chen  and  Miyakoda,  1974;  Jones,  1976;  Madala, 
1977;  and  Miyakoda  and  Rosati,  1977).   They  range  from  in- 
terpolation schemes  to  using  the  finite  difference  equations 
in  providing  the  various  boundary  values.   The  method  used 
in  this  paper  was  the  Dirichlet  condition  with  local  bound- 
ary smoothing  (Chen  and  Miyakoda,  1974). 

The  boundary  points  are  provided  all  necessary  values 
of   U  ,  V  ,  7T  ,  and  T   for  both  one  and  two  way  interaction 
approaches  from  the  solutions  of  the  coarse  mesh  grid. 
These  boundary  values  vary  with  time  as  the  coarse  mesh  ad- 
vances in  time.   Linear  interpolations  in  time  and  space  are 
provided  for  the  fine  grid  where  necessary.   Due  to  the  flux 
form  of  the  advection  terms,  additional  mass  flux  terms 
(ttv  — )   and   (ttu  — )   are  required  outside  the  interface. 
These  values  are  also  interpolated  from  the  coarse  mesh. 

As  described  by  Chen  (1973)  and  Chen  and  Miyakoda  (1974) , 
this  method  of  providing  all  variables  at  the  boundary  is  an 
over specif icat ion.   This  can  cause  a  computational  mode  to 
develop  at  the  boundaries.   A  method  to  reduce  this  computa- 
tional mode,  but  not  to  suppress  significantly  the  physical 
mode  is  to  apply  a  local  smoother  to  all  variables  at  the 
grid  points  next  to  the  boundaries.   Chen  and  Miyakoda  (1974) 
suggested  the  use  of  a  simple  filter  (1-2-1)  applied  normal 
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to  the  boundary  at  the  adjacent  row.   In  this  research  a  five 
point  filter  was  used  to  help  suppress  lateral  waves  along 
the  boundaries.   For  example,  the  smoothing  routine  for  the 
U-component  at  the  west  boundary  is  as  follows: 

U2>J(smoothed)  =  \   U2>J  +  ^^j+D^jffl^j^^.j)  . 

This  filter  was  not  applied  after  every  time  step,  but  after 
every  6-3-3-3  time  steps  (Fig.  3). 

As  with  the  coarse  mesh  boundary  conditions,  the  scheme 
C  grid  arrangement  and  flux  forms  of  the  advective  terms  af- 
fected the  choice  of  the  interface  walls.  A  scheme  C  (I, J) 
grid  point  was  used  for  the  interface  walls.  In  Fig.  7  the 
points  surrounded  by  the  box  are  predictive  points,  whereas 
fine  mesh  lines,  I,  J,  1+9,  J+6  are  interface,  non-predictive 
points . 

Unlike  the  channeled  coarse  mesh  grid,  there  are  no  con- 
servation of  mass  or  energy  requirements  for  the  nested  fine 
mesh.   Therefore,  free  exchange  across  the  interface  is 
necessary.   Providing  the  boundary  values  with  local  smooth- 
ing was  insufficient  to  maintain  free  exchange  of  mass  and 
energy  across  the  interfaces.   This  can  be  seen  in  Fig.  15b, 
which  shows  the  reflections  of  the  U-component  from  the  east 
interface.   A  method  suggested  by  Madala  (1977),  which  con- 
sisted of  calculating  the  mass  convergence  of  the  coarse  mesh 
at  the  interface  to  provide  values  for  certain  boundary  fine 
mesh  mass  flux  terms,  was  attempted  without  success.   (Refer 
to  Appendix  B  for  details  of  the  divergence  boundary  condition.) 
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Problems  arose  at  the  west  and  north  interface  walls  where 
the  horizontal  flux  terms  were  removed  considerably  from 
the  interface.   Changing  the  interface  walls  to  consist  only 
of  the  tt  ,  T  ,  and  $   row  would  have  allowed  better  conser- 
vation of  divergence  between  grids.   This  change,  though, 
would  have  entailed  changing  the  flux  forms  of  the  advective 
terms  and  considerable  indexing  changes.   Time  did  not  permit 
complete  testing  of  this  boundary  condition. 

C.   MARCHING  PROCESS  FOR  THE  NESTED  MODEL 

A  schematic  of  the  Matsuno-leapfrog  time  integration 
scheme  for  both  the  coarse  and  fine  mesh  grids  is  shown  in 
Fig.  3.   Each  sequence  is  45  minutes  in  length.   This  time 
sequence  applies  the  Matsuno  scheme  less  often  in  the  fine 
mesh,  which  results  in  less  damping  of  the  high  frequency 
oscillations.   This  is  in  contrast  to  Chen  and  Miyakoda 
(1974),  who  used  the  Euler-backward  time  step  solely,  and 
Miyakoda  and  Rosati  (1977) ,  who  used  a  time  filter  (Asselin, 
1972)  after  each  time  step. 

Due  to  the  oscillations  generated  at  the  boundaries  of 
the  coarse  and  fine  meshes,  it  was  difficult  to  evaluate  the 
above  time  stepping  scheme.   Computational  modes  did  not 
appear  to  be  significant  in  the  various  model  cases  studied. 
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V.     DESCRIPTION  OF  THE 
EXPERIMENT  AND  RESULTS 

The  main  objective  in  all  nested  grid  models  is  to  show 
that  the  nested  solution  gives  results  similar  to  those  of 
a  uniform  fine  mesh.   This  applies  to  situations  in  which 
the  resolution  of  the  coarse  mesh  grid  is  inadequate  and  the 
fine  mesh  resolution  is  sufficient.   An  attempt  was  made  in 
this  research  to  obtain  a  nested  fine  mesh  solution  for  flow 
over  mountains . 

The  following  is  a  listing  of  the  three  basic  models 
studied  in  this  research. 

(1)  Coarse  mesh  solution: 

a)  without  a  mountain 

b)  with  mountain  A 

c)  with  mountain  B 

(2)  Uniform  fine  mesh  solution: 

a)  without  a  mountain 

b)  with  mountain  A 

(3)  Nested  solution: 

a)  without  a  mountain,  one  way 

b)  without  a  mountain,  two  way 

c)  with  mountain  A,  one  way 

d)  with  mountain  A,  two  way 

Mountain  "A"  is  the  case  in  which  the  mountain  is  built  com- 
pletely across  the  channel,  whereas  mountain  "B"  is  built 
only  half-way  across  the  channel.   Results  to  12  hours  are 
compared,  due  to  problems  developing  after  that  time. 
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A.   COARSE  MESH  SOLUTION 

Hayes  (1977)  studied  the  flow  over  mountains  using  the 
same  primitive  equation  model  (Monaco,  1975),  except  he 
used  the  full  global  model  and  not  a  channeled  version.   Some 
of  his  results  which  can  be  used  in  comparison  are  as  follows: 
(a)  the  mountains  interact  with  the  atmospheric  flow  to  in- 
duce ridging  over  the  mountains  and  troughing  on  the  down- 
stream side;  (b)  damping  of  small  scale  features  with  height; 
and  (c)  increasing  the  grid  resolution  causes  the  scale  of 
the  velocity  component  response  to  decrease,  but  the  magnitude 
remains  about  the  same  when  compared  to  the  coarse  grid 
resolution.   The  grid  size  that  Hayes  used  for  the  coarse 
mesh,  4.5°  long  by  4.0°  latitude,  is  the  same  as  used  by 
this  author. 

Prior  to  testing  the  channeled  model  with  a  mountain,  a 
case  was  computed  without  a  mountain.   This  was  done  to  check 
if  the  initialization  scheme,  which  used  the  linear  balance 
equation,  was  satisfactory  for  the  channeled  model.   The 
model  was  run  to  56  hours  with  the  maximum  variation  in  both 
velocity  components  not  exceeding  +  1  m/sec  over  the  initial 
values.   Some  adjustment  was  expected  because  the  finite 
difference  form  of  the  analytic  solution  is  inexact.   The 
terrain  pressure  does  experience  a  variation  with  time  (Fig. 
3).   This  is  due,  as  explained  in  Chapter  III,  to  the  inte- 
grated mass  flux  in  the  north-south  direction  not  being  equal 
to  zero.   Therefore,  other  than  the  pressure  variation,  the 
model  does  appear  to  be  in  good  balance. 
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Next,  a  meridional  mountain  range  18°  x  750  m  (width  by 
height)  extending  across  the  entire  channel  was  considered. 
Since  the  distance  between  grid  points  in  the  ^-direction  is 
4.5  ,  this  allows  only  three  grid  points  at  each  latitude  to 
represent  the  mountain.   The  middle  channel  U  and  V  compo- 
nents after  six  hours  are  shown  in  Figs.  9A  and  9B  for  three 
levels.   The  mountain  has  developed  to  half  of  its  height  by 
this  time.   Upstream,  the  flow  is  deflected  northward  by  the 
mountain  and  downstream  the  flow  is  deflected  to  the  south. 
At  some  distance  from  the  mountain  both  upstream  and  down- 
stream, variations  due  to  the  mountain  are  less  evident. 
Figure  9B  shows  an  increase  in  the  U-component  in  the  vicinity 
of  the  mountain  top.   The  increase  can  partly  be  explained  by 
the  use  of  the  conservation  of  potential  vorticity  equation, 

g  +  f       ,.  _ 

—. =  constant  , 

Ap 

where  z,      is  relative  vorticity,   f   is  the  Coriolis  para- 
meter and  Ap   is  the  pressure  depth  of  a  column.   In  order 
to  conserve  potential  vorticity  on  the  upstream  side  the  flow 
turns  anticyclonic  and  therefore  an  increase  in  the  V-compo- 
nent  results.   Over  the  mountain  crest  the  flow  turns  cyclonic 
and  the  V-component  decreases,  but  due  to  the  decreased  column 
depth  the  conservation  of  potential  vorticity  is  achieved  by 
an  increase  in  the  U-component. 

The  vertical  velocity  field,  shown  in  Fig.  10,  indicates 
rising  motion  on  the  upstream  side  with  sinking  motion  just 
over  the  mountain  crest  and  downslope.   In  addition,  a 
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secondary  region  of  rising  motion  is  indicated  downstream. 
This  can  be  indicative  of  a  mountain  lee  trough. 

The  temperature  patterns  for  the  middle  of  the  channel 
at  the  six  different  levels  is  shown  in  Fig.  11.   Note  the 
cooling  on  the  upslope  side  of  the  mountain  and  the  warming 
downstream.   This  is  a  direct  result  of  the  vertical  veloc- 
ities induced  by  the  mountain.   It  must  be  noted,  though, 
that  the  temperature  field  over  the  entire  channel  increased 
on  the  average  of  10  K  at  all  levels  from  initial  tempera- 
tures.  This  was  due  in  part  to  term  (1)  of  the  equation  in 
Appendix  A,  the  integrated  n- component  of  momentum  and  the 
imbalance  of  mass  flux  terms  between  the  north  and  south 
boundaries . 

At  12  hours,  the  mountain  has  attained  its  final  height. 
The  horizontal  wind  profiles  for  this  time  are  shown  in 
Figs.  12A  and  B.   The  U-profiles  appear  reasonable  with  the 
highest  zonal  velocity  over  the  mountain  crest,  and  also  a 
decreasing  magnitude  with  height.   The  V-profiles,  though, 
begin  to  show  problems.  At  level  six  the  flow  is  still  being 
deflected  in  the  proper  way,  but  problems  are  appearing  at 
levels  one  and  three.   Although  they,  too,  are  showing 
proper  deflection,  the  maxima  of  level  one  is  the  same  magni- 
tude as  level  six.   This  is  unexpected,  based  on  Hayes1 
results  showing  dampening  with  height. 

The  large  unexpected  maxima  at  level  one  is  possibly  due 

to  two  factors.   The  first  is  due  to  the  method  of  specifying 

4h-  at  the  boundaries  (refer  to  Chapter  III  for  details.) 
Sri 
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The  second  is  the  apparent  development  of  a  meridional,  ver- 
tical oscillation  as  a  result  of  the  vertical  velocity  pat- 
tern induced  by  the  mountains.   Note  in  Fig.  12A  the  apparent 
compensating  northward  flow  downstream  at  level  six  in  re- 
sponse to  the  southward  flow  at  levels  one  and  three.   In 
Fig.  13  the  largest  vertical  velocities  at  level  five,  and, 
therefore,  vertical  transport  of  momentum,  are  shown  to  be 
occurring  along  the  channel  walls  in  the  vicinity  of  the 
mountain.   This  pattern  of  the  largest  vertical  velocities 
occurring  at  the  channel  wall-mountain  intersection  con- 
tinues with  time  and  dominates  the  vertical  velocity  patterns 
after  24  hours. 

To  estimate  the  effect  of  the  mountain- channel  wall  in- 
tersection, a  test  case  was  computed  using  the  half -channel 
mountain  "B".   The  horizontal  velocity  profiles  were  similar 
to  the  channeled  mountain,  except  for  a  slight  increase  in 
the  zonal  flow  north  of  the  mountain.   The  12  hour  vertical 
velocity  pattern,  Fig.  14,  associated  with  level  five  of 
mountain  "B"  is  slightly  different,  especially  along  the 
north  wall,  when  compared  to  the  pattern  of  mountain  "A", 
Fig.  13.   At  12  hours  the  largest  vertical  velocities  appeared 
along  both  the  north  and  the  south  walls.   Even  with  the 
mountain  not  intersecting  the  north  wall,  vertical  velocities 
of  the  same  magnitude  as  in  Fig.  13  developed.   Further  tests 
would  be  required  to  determine  the  minimum  distance  needed 
to  prevent  interaction  of  the  flow  around  the  mountain  and 
the  channel  wall. 
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B.   UNIFORM  FINE  MESH 

The  same  geographic  domain  as  in  the  coarse  mesh  model 
was  used  in  the  fine  mesh  model  with  one- third  the  grid  dis- 
tance.  Due  to  the  increase  in  the  computer  memory  require- 
ments (9X)  and  CPU  (Central  Processing  Unit)  time  (27X) , 
associated  with  this  solution,  few  case  studies  were  com- 
pleted. 

The  case  in  which  there  was  no  mountain  was  computed  to 
evaluate  the  initial  balance  and  pressure,  i.e.  mass,  varia- 
tion over  the  total  field.   Results  through  12  hours  showed 
maximum  velocity  variations  of  +  1  m/sec.   The  average  pres- 
sure over  the  domain  did  vary  with  increases  of  .3  mb  at  six 
hours  and  2.6  mb  at  12  hours.   The  phase  and  period  of  these 
pressure  oscillations  were  difficult  to  calculate  as  the  re- 
quired long  computer  runs  had  to  be  limited  due  to  the  exces- 
sive processing  times.   It  appears,  though,  to  have  a  phase 
close  to  that  of  the  coarse  grids  phase,  but  a  smaller 
magnitude. 

Regrettably,  due  to  time  limitations,  the  uniform  fine 
mesh  with  a  channel  mountain  was  not  run  for  a  sufficiently 
long  time  to  use  in  this  study.   Results  achieved  by  Hayes 
(1977)  in  increasing  the  horizontal  grid  resolution  (£- 
direction  only)  can  be  used  as  a  basis  of  what  to  expect. 
Mainly,  he  found  that  by  increasing  the  resolution  in  the 
^-direction  the  velocity  component  response  decreased,  but 
the  magnitude  of  that  response  remained  approximately  the 
same  when  compared  to  the  coarser  resolution  grid. 
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C.   NESTED  SOLUTION 

Two  methods  of  interaction  are  employed  with  the  nested 
grid  model.   The  first  method  is  one-way  interaction;  values 
from  the  coarse  mesh  are  supplied  to  the  fine  mesh  boundary 
at  appropriate  times,  with  no  corresponding  changes  of  coarse 
mesh  interior  values  from  the  nested  fine  mesh.   The  other 
method  is  two-way  interaction.   Although  boundary  values  for 
the  fine  mesh  are  still  provided  by  the  coarse  mesh,  replace- 
ment of  coarse  mesh  values  by  interior  nested  fine  mesh 
values  at  coincident  points  is  done  at  appropriate  times. 
For  example:   v(I+3,J+2)  of  fine  mesh  replaces  v(I+2,  J+2) 
of  coarse  mesh.   A  five  point  averaging  scheme  surrounding 
each  fine  mesh  coincident  point  is  used  to  calculate  the 
replacement  fine  mesh  value. 

Prior  to  testing  the  nested  grid  with  a  mountain  present, 
one  case  with  each  type  of  interaction  was  run  without  a 
mountain  to  check  the  initial  balance.   As  mentioned  in  sec- 
tion 1  of  this  chapter,  the  coarse  mesh  is  in  good  balance. 
The  v  and  u  components  in  the  middle  region  of  the  fine  mesh 
after  six  hours  are  shown  in  Figs.  15A  and  B.   The  corre- 
sponding coarse  mesh  solution  can  be  seen  at  either  end  of 
each  plot.   The  v-component  profile  shows  level  l(a  =  .08) 
to  have  the  largest  value.   This  can  be  due  to  the  following: 

1)  the  oscillations  of  the  height  surfaces  in  the  north- 
south  direction  forced  by  the  coarse  mesh  boundary  values, 

2)  the  interface  boundary  conditions,  and  3)  the  phase  dif- 
ferences in  pressure  oscillations.   The  u-component  profile, 
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Fig.  15B,  shows  the  reflection  of  energy  from  the  east  wall. 
Insufficient  interface  boundary  conditions  did  not  allow 
free  exchange  of  mass  and  energy  across  the  interface, 
thereby  resulting  in  the  reflections. 

In  addition,  Fig.  5  shows  the  pressure  variations  for 
the  nested  (one-way)  fine  mesh.   These  oscillations  appear 
to  follow  the  variations  of  the  1.5°  uniform  fine  mesh 
(nested  domain)  spacing  out  to  eight  hours.   After  that  time, 
the  coarse  mesh  variations  force  the  nested  variations .  This 
could  lead  to  imbalances. 

When  the  two-way  interaction  scheme  is  used,  the  fine 
mesh  has  significant  influence  on  the  coarse  mesh,  Fig.  16. 
In  one-way  interaction  the  coarse  mesh  v-profile  would  show 
only  a  +.5  m/sec  value  for  the  three  levels.   It  should  be 
noted,  though,  that  the  variations  of  the  fine  mesh  hori- 
zontal velocity  components  were  not  as  large  in  the  two-way 
as  in  the  one-way.   This  was  due  to  the  transfer  of  energy 
out  of  the  nested  grid  into  the  coarse  grid.   This  allowed 
any  oscillations  present  in  the  fine  mesh  region  to  propa- 
gate to  the  coarse  mesh  region.   These  oscillations,  though, 
aggravated  the  problems  at  the  coarse  mesh  boundary  wall. 
From  12  hours  and  beyond,  the  oscillations  in  both  grids 
become  unrealistically  large. 

The  last  series  of  experiments  was  computed  using  the 
nested,  two-way  interacting  model  with  a  18  x  750  m  mountain 
present.   This  allows  three  grid  points  in  the  coarse  mesh 
and  11  grid  points  in  the  fine  mesh  to  represent  the  mountain 
at  each  latitude.   The  resulting  horizontal  velocity  profiles 
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after  six  hours  for  the  fine  and  coarse  meshes  are  shown  in 
Figs.  17  and  18.   The  mountain  has  developed  to  half  its 
height  by  this  time.   The  v-profiles  in  Fig.  17A  show  the 
general  upstream  northward  deflection  and  downstream  south- 
ward deflection.   In  addition,  the  upstream  flow  appears  to 
be  deflected  northward  only  in  the  immediate  vicinity  of  the 
mountain.   This  is  different  than  the  coarse  mesh  solution 
which  shows  deflections  three  gird  points  upstream  from  the 
mountain.   This  narrowing  of  the  velocity  component  response 
compares  favorably  with  that  of  Hayes  (1977) .   One  result 
which  does  not  compare  is  the  lack  of  dampening  with  height. 
This  is  probably  a  result  of  the  improper  v-component  values 
from  the  coarse  mesh  being  forced  at  the  nested  boundary  in- 
terface. 

The  u-profile  (Fig.  17B) ,  for  level  6(a  =  .92)  does  indi- 
cate an  increase  in  the  zonal  flow  over  the  mountain  top,  and 
then  a  return  to  the  upstream  value.   The  upper  two  levels 
(a  =  .42,  .08),  though,  show  an  increase  in  the  zonal  flow 
downstream  of  the  mountain.   This  is  probably  caused  by  the 
reflections  from  the  east  interface  wall. 

The  coarse  mesh  velocity  profiles  are  shown  in  Figs.  18A 
and  B.   When  compared  to  Figs.  9A  and  B,  the  one-way  inter- 
acting mountain  case,  few  significant  differences  are  noticed. 
The  v-profiles  show  the  proper  deflection  north  and  southward, 
with  the  zonal  flow  increasing  in  the  region  of  the  mountain 
top.   Since  the  flow  in  the  region  of  the  mountains  is  just 
an  average  of  the  nested  region  flow,  the  fine  mesh  solution 
must  have  been  good  for  the  first  six  hours. 
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In  Figs.  19  and  20  are  shown  the  v-coinponent  profiles 
and  vertical  velocities  for  the  coarse  mesh  at  12  hours. 
The  v-profiles  show  a  large  level  l(a  =  .08)  northward  com- 
ponent with  a  compensating  lower  level  G(a  =  .92)  southward 
component.   The  fine  mesh  v-profile  is  similar,  except  the 
deflections  are  spread  wider  over  its  domain.   Part  of  the 
reason  for  the  large  upper  level  v-component  is  the  inter- 
action of  the  mountain  and  the  channel  walls.   As  mentioned 
in  section  1,  large  vertical  velocities  are  created  at  the 
intersection  with  channel  wall  points,  which  causes  an 
artificial  meridionally  oriented,  vertical  circulation.   The 
same  circulation  apparently  occurs  in  the  fine  mesh,  due  to 
interface,  mountain  intersection.   This  energy  is  then  par- 
tially transferred  to  the  coarse  mesh  via  the  two-way  inter- 
action.  The  coarse  mesh  vertical  velocities  are  shown  in 
Fig.  20.   The  fine  mesh  influence  is  immediately  apparent. 
By  this  time,  the  solutions  to  the  nested  grid  model  are 
dominated  by  unrealistic  oscillations.   Similar  results  were 
obtained  in  the  one-way  interaction  solutions. 
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VI.   SUMMARY  AND  CONCLUSIONS 

This  thesis  was  undertaken  to  study  the  flow  over  moun- 
tains using  a  channeled  model  with  an  interior  nested  grid. 
The  work  involved  two  stages.   The  first  was  to  develop  an 
east-west  channel  model  from  a  global  model.   The  second  was 
the  nesting  of  a  finer  resolution  grid. 

In  channeling  the  modified  Arakawa  global  model  by 
Monaco  and  Williams  (1975) ,  oscillations  and  noise  arose 
along  the  channel  walls.   The  scheme  C  grid  arrangement  and 
flux  form  of  the  advective  terms  lead  to  certain  boundary 
conditions.   These  proved  to  be  adequate  when  no  disturbance 
was  present,  but  insufficient  when  either  a  mountain  and/or 
nested  grid  were  introduced.   Difficulties  also  arose  at  the 
mountain- channel  wall  intersection. 

In  nesting  a  finer  resolution  grid  several  things  were 
learned.   The  numerical  response  in  the  ^-direction  decreases 
with  the  finer  resolution  when  compared  to  the  coarse  mesh 
solution.   In  addition,  the  two-way  interaction  does  allow 
partial  transfer  of  energy  to  the  coarse  mesh.   Several  prob- 
lems were  also  encountered. 

The  grid  arrangement  and  flux  form  of  advective  terms 
lead  to  similar  boundary/ interface  conditions  as  in  the 
coarse  mesh.   Difficulties  did  not  arise  at  the  interface 
walls  due  to  extrapolation  methods  in  providing  boundary 
values,  but  they  resulted  due  to  not  allowing  free  exchange 
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of  energy  and  mass  across  the  interface  walls.   This  caused 
internal  reflected  fine  mesh  waves,  which  eventually  de- 
stroyed the  fine  mesh  solution.   The  mountain- interface 
wall  intersection  also  accounted  for  some  of  the  problems. 
As  in  the  coarse  mesh,  large  vertical  velocities  developed 
at  these  intersection  points,  which  are  partly  responsible 
for  the  vertical  meridional  oscillations. 

Based  on  the  above  results,  the  following  is  recommended: 

(a)  expand  the  channel  in  the  north-south  direction  suffi- 
ciently to  allow  isolating  the  mountain.   This  would 
avoid  the  problems  caused  by  the  mountain  channel  wall 
intersection; 

(b)  simplify  the  flux  forms  in  the  primitive  equations  for 
use  over  the  whole  domain; 

(c)  modify  the  channel  walls  to  consist  of  one  latitude 
$  ,  t  ,  u  row.   This  leads  to  an  unequal  predictive 
region  for  the  v-component; 

(d)  modify  the  interface  walls  in  the  nested  grid  to  consist 
of  one  latitude  (north  and  south  walls)  and  one  longi- 
tude (east  and  west  walls)  row.   This  leads  to  unequal 
predictive  rows  for  v  and  u- components ,  respectively; 

(e)  alter  linear  interpolation  scheme  to  higher  order  scheme 
for  providing  the  nested  grid  boundary  with  coarse  mesh 
values . 

Although  this  thesis  was  not  entirely  successful  in  mathe- 
matically simulating  flow  over  mountains,  it  did  identify  and 
detail  some  specific  problems.   It  is  hoped  that  this  study 
will  prove  useful  to  those  attempting  a  similar  study  in  the 

future . 
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APPENDIX  A 
NORTH  AND  SOUTH  BOUNDARY  CONSIDERATIONS 

The  constraint  of  no  net  mass  flux  in  the  north-south 
direction  is  a  primary  factor  to  consider  when  determining 
boundary  conditions  for  the  north  and  south  walls.   Due  to 
the  cyclic  condition  in  the  east-west  direction,  the  only 
equation  of  motion  needed  to  be  considered  is  equation  (2-3) , 
the  v- component  equation: 

Tt(Sr  V)  +  TC(1T  v)  +  ^("m-  v)  +  To<mi   v) 
(1)       (2)        (3)        (4) 

flTU  ,31  3         lx  7Tr9<I>  „         dTT .  TT 

+  ^T+(vhH-u^T  m")7TU  +  E[37T  +  aa    87?1    =  mn"  Fn    ' 
(5)  (6)  (7)  (8) 

To  assure  no  net  mass  flux,  the  above  equation,  integrated 
in   £   and   a  ,  must  balance. 

Following  is  an  analysis  of  each  term,  and  what  must  be 
done  to  accomplish  the  constraint  condition. 
Term  (1) : 


//  &&**>* -h   // 


(— )dad£  =  0 


TTV 

In  the  case  where  no  mountain  is  present,   —  starts  at 
zero  and  should  remain  zero  as  long  as  the  model  initiali- 
zation is  correct.   When  a  mountain  is  present,  there  must 
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be  a  shift  of  mass  to  another  region  in  order  to  satisfy 
the  integral. 
Term  (2) : 


J    ir(ir  v) 


dadC  =  0 

This   is   due    to   the   cyclic   condition  east-west. 
Term    (3) : 

|_(H  v)dodc  .    / /v^do*  +  //=£  dadC 

The  first  term  on  the  right  hand  side  is  taken  care  of  by 
reflecting  the  normal  mass  flux  quantities  along  the  bound- 
aries.  In  the  no  mountain  case,  the  second  term  starts  at 
zero  and  remains  close  to  zero  as  the  model  advances  in  time. 
When  a  mountain  is  introduced,  this  term  becomes  non-zero. 
This  problem  is  further  accentuated  by  the  imbalance  in  the 
pressure  gradient  term  (7) . 
Term  (4) : 


// 


3  (I£  v)dad?  =  0 


If  a  =  0   at   a  =  1  and   a  =  0  ,  this  term  is  satisfied, 
Term  (5) : 


// 


Ettu  dad?   +    0 


mn 
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This  term  is  not  equal  to  zero,  but  will  primarily  be  bal- 
anced by  the  pressure  gradient  term.   The  term  will  further 
vary  if  a  mountain  is  present. 
Term  (6) : 


J  }<rkk-»hh™**-   jf 


7  |g.  I  dad? 


S     a  5o 


-// 


Alr5dod5  +  0 


The  first  term  on  the  right  hand  side  starts  at  zero  (v=0) , 
and  will  remain  close  to  zero  as  the  model  advances  in  time. 
The  second  term  is  not  zero,  because  the  integrated  scaling 
factor  is  non-zero.   The  contribution,  though,  is  very  small 
compared  to  the  other  terms  in  the  equation.   Its  importance, 
however,  increases  as  the  velocity  components  increase. 
Since  the  other  terms  are  not  exactly  zero,  this  term  adds 
to  the  imbalance. 
Term  (7) : 


f  I       -[^ 
J   J   ml8n 


+  act  p-]doa^   *  0 


This  term  is  non-zero  and  primarily  balanced  by  the  £-compo- 
nent  of  velocity.   In  order  to  prevent  mass  variations, 
though,  its  initial  balance  along  the  boundaries  must  be 
maintained.   The  initial  terrain  pressure  gradient  was  rela- 
tively simple  to  maintain  but  the  initial  geopotential  height 

gradient  was  not  so  simple. 
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Term    (8)  : 

dodE,   =   0 


J       J       mn      T) 


5        a 


Tliis  model   is   frictionless 
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APPENDIX  B 
DIVERGENCE  BOUNDARY  CONDITION 

Free  exchange  of  mass  and  energy  across  the  nested  grid 
interface  is  a  necessity  to  prevent  reflections  from  the 
walls  destroying  the  fine  mesh  solutions.   In  addition  to 
the  two-way  interaction  scheme,  a  divergence  boundary  condi- 
tion after  Iladala  (1977)  was  attempted. 

This  boundary  condition  entailed  calculating  the  coarse 
mesh  divergence  along  the  interface,  and  interpolating  spa- 
tially and  temporally  to  determine  the  fine  mesh  mass  flux 
quantity  on  the  outside  of  the  interface.   The  following 
formulas  were  used  for  the  east  and  south  interfaces  respec- 
tivelv: 


<™  fh,J  -  DIVcoarse  +  <™  §>I-1,J  "  <™  #>I, Jfl 


+  (ttv^)t    _ 
v        nrI,J 


to  £>I,J  =  '^coarse  +  <™  f  \j  -  (™  f  >I.1,J 

The  interface  divergence  boundary  conditions  were  not  applied 
along  the  north  and  west  interfaces,  due  to  the  problems  men- 
tioned in  chapter  IV.   With  this  boundary  condition  used 
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along  the  east  and  south  walls,  results  were  still  unsatis- 
factory.  Reflections  were  still  occurring,  indicating  a 
free  transfer  of  mass  and  energy  was  not  present. 

Madala  (1977)  apparently  solved  this  problem  by  using 
weighted  extrapolated  fine  mesh  U  and  V  values,  along  with 
coarse  mesh  values,  in  determining  the  mass  flux  divergence 
terms  for  the  fine  mesh  boundary.   This  modified  the  forc- 
ing of  the  coarse  mesh  divergence  fields  at  the  interface. 
It  should  be  noted  that  Madala  used  only  the  tt  ,  T  ,  and 
(p      row  of  the  scheme  C  grid  for  the  boundary.   In  addition, 
he  used  a  simplified  flux  form,  thereby  he  did  not  need  to 
specify  as  many  external  boundary  points.   This  makes  the 
north  and  west  interface  conditions  exactly  the  same  as  the 
east  and  south,  whereas  in  this  research  they  would  have 
been  different. 
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nested  grid  smoother 


Figure  8.   Time  step  routines  for  a  45  minute  sequence. 
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Figure  9B. 


U-components  for  coarse  mesh  with  mountain  A. 
time  =  6  hours 
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U-components  for  coarse  mesh  with  mountain  A, 
time  =  12  hours 

a  -    .08  (solid  line)  a  =    .42  (dashed  line) 
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Figure  15B. 


U-components  for  fine  mesh  (one-way)  with  no 
mountain. 

time  =  6  hours 
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V- components  for  coarse  mesh  (two-way)  with 
no  mountain. 

time  =  6  hours 
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a  =  .92  (dotted  line) 
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Figure  17B.   U-components  for  fine  mesh  (two-way)  with 
mountain  A. 

time  =  6  hours 
a  =  o08  (solid  line)   a  =  . 42  (dashed  line) 

a   =  .92  (dotted  line) 


65 


10r 


V  m 
sec 


-5 


-ry^   I 


-10*- 

Figure  18A, 


20r 


V- components  for  coarse  mesh  (two-way)  with 
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Figure  18B.   U-components  for  coarse  mesh  (two-way)  with 
mountain  A. 

time  =  6  hours 
a  =  .08  (solid  line)   a  =  .42  (dashed  line) 
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